
#include "utility.h"
#include <Eigen/Dense>
#include <cassert>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
constexpr double re_num{100000.0};
void OutputResult(const VarArray<Variable, 2, 0> &flow,
                  const VarArray<Point, 2, 0> &grid, ostream &os);
int main() {

  VarArray<Point, 2, 0> grid;
  string grid_file_name{"flat_plate2.x"};
  ReadMeshASSCII(grid_file_name, grid);

  int npoint_in_ydir = grid.Dim(1);
  int npoint_in_xdir = grid.Dim(0);
  //   cout << npoint_in_ydir << endl;
  // 注意
  VarArray<Variable, 2, 0> flow(npoint_in_xdir, npoint_in_ydir);
  InitializeFlow(flow);

  CalcV update_v(flow, grid);
  update_v.CalcCoefficient();
  //  迭代求解
  for (int i = 1; i < npoint_in_xdir; ++i) {

    InitializeNextLine(flow, i);

    int count = 0;
    double tmpu, tmpv, res_u, res_v, res_u1, res_v1;

    do {

      count++;
      // // 计算之前未知时间层需要使用已知时间层初始化
      res_u = NewtonIterForU(flow, i, grid, re_num);
      update_v.CalcFxk(i);
      // res_v = update_v.UpdateVFourthOrder(i);
      //   res_v = UpdateVfirstOrder(flow, i, grid);
      // res_v = UpdateVthirdOrder(flow, i, grid);
      //   res_v = UpdateVsecondOrder2(flow, i, grid);

      if (count == 1) {
        res_u1 = res_u + 1.0e-28;
        res_v1 = res_v + 1.0e-28;
      }

      tmpu = log10(res_u / res_u1);
      tmpv = log10(res_v / res_v1);

    } while (tmpu > -12.0 && tmpv > -12.0 && res_u > 1.0e-14 &&
             res_v > 1.0e-14 && count < 100);

    cout << "counti = " << count << endl;
    cout << "tmpu = " << tmpu << endl;
    cout << "tmpv = " << tmpv << endl;
    // cout << "i = " << i << endl;
    // exit(EXIT_FAILURE);
  }

  cout << scientific << setprecision(14);
  for (int j = 0; j < npoint_in_ydir; ++j)
    cout << flow(npoint_in_xdir - 1, j).u << endl;
  cout << "=============================" << endl;
  for (int j = 0; j < npoint_in_ydir; ++j)
    cout << flow(npoint_in_xdir - 1, j).v << endl;

  //   ofstream fluid{"flow.plt"};
  //   OutputResult(flow, grid, fluid);
  //   fluid.close();
  return 0;
}

void OutputResult(const VarArray<Variable, 2, 0> &flow,
                  const VarArray<Point, 2, 0> &grid, ostream &os) {
  int npoint_in_ydir = flow.Dim(1);
  int npoint_in_xdir = flow.Dim(0);

  os << "TITLE = \" visualization of the 2Dflow\" \n"
     << "VARIABLES = \"x\", \"y\", \"u\", \"v\"\n";

  os << "ZONE  I=" << npoint_in_xdir << ", J=" << npoint_in_ydir
     << ", DATAPACKING=POINT\n";

  os << std::scientific << std::setprecision(16);

  for (int j = 0; j < npoint_in_ydir; ++j) {
    for (int i = 0; i < npoint_in_xdir; ++i) {
      os << grid(i, j).x << " " << grid(i, j).y << " " << flow(i, j).u << " "
         << flow(i, j).v << "\n";
    }
  }

  // string last{"myBlasius.plt"};
  ofstream blasius{"myBlasius.plt"};
  blasius << std::scientific << std::setprecision(16);
  blasius << "TITLE = \" visualization of the last line\" \n"
          << "VARIABLES = \"y\", \"u\", \"v\", \"eta\"\n";

  blasius << "ZONE  I=" << npoint_in_ydir << endl;

  for (int j = 0; j < npoint_in_ydir; ++j) {
    const Point &p = grid(npoint_in_xdir - 1, j);
    blasius << p.y << " " << flow(npoint_in_xdir - 1, j).u << " "
            << flow(npoint_in_xdir - 1, j).v << " "
            << p.y / p.x * sqrt(re_num * p.x) << endl;
  }
  blasius.close();
}